---
title: "Replication for 'How to Stop Contagion: Applying Network Science to Evaluate the Effectiveness of Covid-19 Vaccine Distribution Plans'"
header-includes:
author: 
     - Olga V. Chyzh^[Assistant Professor, University of Toronto, olga.chyzh@utoronto.ca] 
date: '`r format(Sys.time(), "%B %d, %Y")`'
output: html_document
---



```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, message=FALSE, warning=FALSE, eval=TRUE,fig.pos = "!t", out.extra = "")

library(kableExtra)
library(sna)
library(igraph)
library(ggplot2)
library(gridExtra)
library(ggpubr)
library(tidyverse)
library(magrittr)
library(ergm)
library(lubridate)
library(stringr)
library(ggplot2)
library(zoo)
library(cem)
library(optmatch)
library(RItools)
library(survival)
library(ergm)
library(dplyr)
library(tidyr)
library(broom)

library(modelsummary)
source("./scripts/plot_round2.R")
source("./scripts/run_combine_tabs.R")
source("./scripts/sum_net_degree.R")
```


## Table 1: A Comparison of States with Early/Late Vaccine Prioritization for Grocery Workers

```{r}
d<-read.csv("./data/state_data1.csv", header=TRUE) %>% mutate(grocery_elig=ymd(grocery_elig))
mybal<-xBalance(gro_priority ~ Gov_Rep+log(GDP_cap)+median_age1519+log(totpop1519)+I(medinc1519/1000)+perc_unemploy1519+percBA1519+urb2010+percblack1519+perclatino1519+biden_marg+log(Cum_covid),data=d, 
  report = c("adj.means","p.values","chisquare.test"))
mybal
```


## Figure 1: Centrality in the Interaction Network

```{r degree-fig, fig.cap = "Centrality in the Interaction Network", fig.align="center", out.width="50%", echo=TRUE}
set.seed(1708)
data(coleman)
ggraph<-graph_from_adjacency_matrix(coleman[2,,], mode="undirected", diag=FALSE)
ggraph <- igraph::delete.vertices(ggraph , which(igraph::degree(ggraph)==0))
degree.cent <- igraph::degree(ggraph)
ggplot() + geom_histogram(aes(x=degree.cent), fill="grey", colour="black")+xlab("Degree")+ylab("Count")+theme_bw()
```

## Figure 2: Contagion in the Interaction Network

```{r contagion-fig, fig.cap = "Contagion in the Interaction Network", fig.align="center", out.width="50%", echo=TRUE}

mynet<-(as.matrix(igraph::as_adj(ggraph)))
patient0d<-names(sort(degree.cent, decreasing=TRUE))[1]
plot_round2(mynet=mynet,patient0=patient0d,myname=" ")
```


## Figure 3: Vaccination Scenarios

```{r vaccine-fig1, fig.cap = "Vaccination Scenarios", out.width="80%", fig.align="center", echo=TRUE}
par(mfrow=c(1,2),mar = c(.5, .5, 1.5, .5)) 
mynet<-(as.matrix(igraph::as_adj(ggraph)))
patient0<-names(sort(degree.cent, decreasing=TRUE))[1]

n=10
vaccine_low<-as.data.frame(degree.cent) %>% filter(row.names(.)!=patient0) %>% slice_min( degree.cent,n=n, with_ties=FALSE) %>% row.names()
vaccine_centr<- as.data.frame(degree.cent) %>% filter(row.names(.)!=patient0) %>% slice_max( degree.cent,n=n, with_ties=FALSE) %>% row.names()

vulnerable<-vaccine_low


plot_round2(mynet=mynet,patient0=patient0,vaccinated=vaccine_low,vulnerable=vaccine_low,myname="Low Degree")
plot_round2(mynet=mynet,patient0=patient0,vaccinated=vaccine_centr,vulnerable=vaccine_low,myname="High Degree")


```


## Table 2: Round 2 Summary for 10,000 Simulations of the Interaction Network among 73 Individuals

```{r sim-res, echo=TRUE, eval=TRUE}
friends<-network(mynet, matrix.type="adjacency", 
                directed=FALSE, loops=FALSE)
m1<-readRDS("./data/ergm1.rds")
set.seed(1314)
g.sim <- simulate(m1, nsim=10000, 
                  basis=friends, control=control.simulate.ergm(
                   MCMC.burnin=1000,
                   MCMC.interval=100))

myres<-lapply(1:length(g.sim), function(i){
  sum_net_degree(mynet=g.sim[[i]]) })
arr <- array( unlist(myres) , c(3,4,length(g.sim)) )
saveRDS(arr,"./data/sim_res.rds")

mytable<-as.data.frame(as.matrix(NA,3,4))
for (i in 1:3){
  for (j in 1:4){
   m<-format(round(mean(arr[i,j,]),2),nsmall=2)
   se<-format(round(sd(arr[i,j,]),2),nsmall=2)
   mytable[i,j]<-paste0(m, " (",se,")")
  }
}

mytable[,2]<-c(0,10,10)
mytable[2,4]<-0


row.names(mytable)<-row.names(myres[[1]])
names(mytable)<-names(myres[[1]])
saveRDS(mytable,"./data/mytable.rds")

mytable<-readRDS("./data/mytable.rds")
mytable<-mytable[c(2,1,3,4)]
mytable %>% kable(caption = "Round 2 Summary for 10,000 Simulations of the Interaction Network among 73 Individuals",align='cccc') %>%
  kable_paper("hover", full_width = F) %>%
   kableExtra::footnote(general="Cell values are means over 10,000 simulations. Numbers in parentheses are standard deviations. \\\\textit{Not Infected} does not include \\\\textit{Vaccinated}.", escape=FALSE,threeparttable=TRUE)
```







## Table 3: Matched Counties

```{r }
mydata<-read.csv("./data/CAOR_data.csv", header=TRUE) %>% mutate(Date=ymd(Date)) %>% filter(Date<ymd("2021-03-29")+weeks(2) & Date>ymd("2020-12-16"))

mydata %>% filter(Date==ymd("2020-12-17")) %>% 
  select(State, County, lnCumCovid1000_pre=lnCumCovid1000,CumCovid_pre=Cum_covid) %>%  
  left_join(mydata, by=c("State","County")) -> mydata


tomatch<-mydata %>%  filter(Date==ymd("2020-12-17")) %>% select(State, County, GDP_co_2019,ln_GDP_co_2019,ln_totpop1519,totpop1519,perc_unemploy1519,percBA1519,urb2010,percblack1519,perc_otherrace1519,perclatino1519,perc_foreign1519,Bidenmarg, grocery_wkrs100, age_65plus_prop, indoor_dine,
                           CumCovid_pre,lnCumCovid1000_pre) 
todrop<-c("State","County","totpop1519","GDP_co_2019","CumCovid_pre","grocery_wkrs100", "age_65plus_prop", "indoor_dine")
mat1 <- cem(treatment = "State",
            cutpoints=list(ln_GDP_co_2019=4,
                           ln_totpop1519=6,
                           perc_unemploy1519=2,
                           percBA1519=2, 
                           urb2010=2,
                           percblack1519=3,
                           perc_otherrace1519=3,
                           perclatino1519=2, 
                           perc_foreign1519=2,
                          Bidenmarg=4),
            baseline.group = "Oregon",
            drop=todrop,data = tomatch, keep.all = TRUE)

matches<-tomatch[mat1$matched==TRUE,] %>% mutate(matched="TRUE")
matched_data1<-mydata %>% left_join(matches, by=c("State","County")) %>% filter(matched==TRUE)

CA<- tomatch %>% ungroup() %>% 
  mutate(strata=mat1$strata, matched=mat1$matched) %>% 
  filter(matched=="TRUE") %>% dplyr::select(State, County, strata) %>% filter(State=="California")

matched_counties <- tomatch %>% ungroup() %>% 
  mutate(strata=mat1$strata, matched=mat1$matched) %>% 
  filter(matched=="TRUE") %>% dplyr::select(State, County, strata) %>% filter(State=="Oregon") %>% left_join(CA, "strata")

matched_counties

```

## Table 4: Balance Between Oregon and California in the Matched Sample

```{r}
mybal<-xBalance(as.numeric(as.factor(State))~ln_GDP_co_2019+ln_totpop1519+perc_unemploy1519+percBA1519+urb2010+percblack1519+perclatino1519+perc_otherrace1519+perc_foreign1519+Bidenmarg+age_65plus_prop+ indoor_dine+lnCumCovid1000_pre,data=matches, 
  report = c("chisquare.test","std.diffs","p.values"))

mybal$results
mybal$overall

```

## Table 5: The Effect of Vaccine Eligibility to Grocery Employees on New Daily Covid-19 Cases (logged)

```{r}
matches<-tomatch[mat1$matched==TRUE,] %>% mutate(matched="TRUE") %>% select(State,County,matched)


matched_data1<-mydata %>% left_join(matches, by=c("State","County")) %>% filter(matched==TRUE)

mydata$State_grocery_day<-mydata$State_dum*mydata$grocery_day
matched_data1$State_grocery_day<-matched_data1$State_dum*matched_data1$grocery_day
m1_all<-lm(ln_Newcases7~grocery_day+State_dum+State_grocery_day+lnCumCovid1000_pre+ln_GDP_co_2019+ln_totpop1519+perc_unemploy1519+percBA1519+urb2010+percblack1519+perclatino1519+perc_otherrace1519+perc_foreign1519+Bidenmarg+ age_65plus_prop+ indoor_dine+lag3_lnNewCasescap1000, data=mydata)

m1_matched<-lm(ln_Newcases7~grocery_day+State_dum+State_grocery_day+lnCumCovid1000_pre+ln_GDP_co_2019+ln_totpop1519+perc_unemploy1519+percBA1519+urb2010+percblack1519+perclatino1519+perc_otherrace1519+perc_foreign1519+Bidenmarg+ age_65plus_prop+ indoor_dine+lag3_lnNewCasescap1000, data=matched_data1)
#options("modelsummary_format_numeric_latex" = "plain")

modelsummary(list("Full Sample"=m1_all, "Matched Sample"=m1_matched), estimate="{estimate}{stars} ({std.error})", statistic=NULL, 
             gof_omit='logLik|BIC|AIC',
          fmt=3,
          caption="The Effect of Vaccine Eligibility to Grocery Employees on New Daily Covid-19 Cases (logged)",
          coef_rename=c("grocery_day"="Day of Treatment", 
                       "State_dum"="California",
                       "State_grocery_day"="California*Day of Treatment",
                       "lnCumCovid1000_pre"="Cumulative Cases, logged",
                       "ln_GDP_co_2019"="County GDP, logged",
                       "ln_totpop1519"="County Population, logged",
                       "perc_unemploy1519"="Unemployment Rate",
                       "percBA1519"="Percent BA Degree",
                       "urb2010"="Urbanization",
                       "percblack1519"="Percent Black",
                       "perclatino1519"="Percent Latino",
                       "perc_otherrace1519"="Percent Other Race",
                       "perc_foreign1519"="Percent Foreign",
                       "Bidenmarg"="Biden's Margin",
                       "grocery_wkrs100"="Grocery Wkrs, per 100 res.",
                       "age_65plus_prop"="Prop. Aged 65+",
                       "indoor_dine"="Indoor Dining Ban",
                       "lag3_lnNewCasescap1000"="Lagged DV",
                       "(Intercept)"="Constant"))
```


## Figure 4: Marginal Effect of Vaccine Eligibility to Grocery Workers. 

```{r fig-marg, fig.cap="Marginal Effect of Vaccine Eligibility to Grocery Workers. Error bars represent 95\\% CIs."}

grocery_day=seq(from=min(mydata$grocery_day),to=max(mydata$grocery_day), by=1)

marg_state<-m1_all$coef[3]+m1_all$coef[4]*grocery_day
se_state<-sqrt(vcov(m1_all)[3,3]+vcov(m1_all)[4,4]*(grocery_day^2)+2*grocery_day*vcov(m1_all)[3,4])

p<-ggplot(data=cbind.data.frame(grocery_day,marg_state,se_state)) + geom_point(aes(x=grocery_day, y=marg_state))  + 
  geom_errorbar(aes(x=grocery_day, ymin=marg_state-2*se_state, ymax=marg_state+2*se_state), width=.1) +
  geom_hline(aes(yintercept=0))+
  scale_y_continuous("Marginal Effect of California", limits=c(-.8,.3)) +
  ggtitle("Full Sample (Model 1)")+
  scale_x_continuous("Day of Treatment", limits=c(0,20))+
  theme_bw()+
  theme(plot.title = element_text(size = 9, face = "plain"))

# (Matched Sample):
grocery_day=seq(from=min(matched_data1$grocery_day),to=max(matched_data1$grocery_day), by=1)

marg_state<-m1_matched$coef[3]+m1_matched$coef[4]*grocery_day
se_state<-sqrt(vcov(m1_matched)[3,3]+vcov(m1_matched)[4,4]*(grocery_day^2)+2*grocery_day*vcov(m1_matched)[3,4])

p1<-ggplot(data=cbind.data.frame(grocery_day,marg_state,se_state)) + geom_point(aes(x=grocery_day, y=marg_state))  + 
  geom_errorbar(aes(x=grocery_day, ymin=marg_state-2*se_state, ymax=marg_state+2*se_state), width=.1) +
  geom_hline(aes(yintercept=0))+
  scale_y_continuous("", limits=c(-.8,.3)) +
  ggtitle("Matched Sample (Model 2)")+
  scale_x_continuous("Day of Treatment", limits=c(0,20))+
  theme_bw()+
  theme(plot.title = element_text(size = 9, face = "plain"))

myfig<-ggarrange(p+ rremove("ylab") + rremove("xlab"),
                 p1+ rremove("ylab") + rremove("xlab"),
                 ncol=2, nrow=1, common.legend = TRUE, legend="bottom")
annotate_figure(myfig, left = "Marginal Effect of California",
                bottom = "Day of Treatment")

```

## Figure 5: Temporal Trends (State Averages) in Covid-19 Cases, December, 2020–April, 2021.

```{r CAOR-cases, eval=TRUE, fig.cap="Temporal Trends (State Averages) in Covid-19 Cases, December, 2020--April, 2021.", results = FALSE}
p<-mydata %>% 
  group_by(State, Date) %>% 
  summarise(NewCases1000=100000*sum(Newcases7)/sum(exp(ln_totpop1519)),.groups="keep") %>%
  ungroup() %>%
  ggplot(aes(x=Date,y=NewCases1000, group=State, color=State, linetype=State))+
  geom_line(linewidth=1)+
  geom_vline(aes(xintercept=ymd("2021-03-01")), linetype=1, color="gray20")+
  geom_vline(aes(xintercept=ymd("2021-03-01")+weeks(2)), linetype=2, color="gray20") +
  geom_text(aes(x=ymd("2021-03-01"), label="grocery wkrs eligible in CA", y=35), colour="gray50", angle=90, vjust = -.5, size=3)+
  geom_text(aes(x=ymd("2021-03-01")+weeks(2), label="14 days after grocery wkrs eligible in CA", y=35), colour="gray50", angle=90, vjust = -.5, size=3)+
  scale_y_continuous("New Covid-19 Cases per 100,000 residents, logged scale", trans="log2", limits=c(5,120))+ xlab("Full Sample")+
  scale_color_manual(values=c("cornflowerblue","red"))+
  theme_bw()+
  theme(axis.text=element_text(size=10),
        axis.title=element_text(size=10))


p1<-matched_data1  %>% 
  group_by(State, Date) %>% 
  summarise(NewCases1000=100000*sum(Newcases7)/sum(exp(ln_totpop1519)),.groups="keep")%>%
  ungroup() %>%
  ggplot(aes(x=Date,y=NewCases1000, group=State, color=State, linetype=State))+
  geom_line(linewidth=1)+
  geom_vline(aes(xintercept=ymd("2021-03-01")), linetype=1, color="gray20")+
  geom_vline(aes(xintercept=ymd("2021-03-01")+weeks(2)), linetype=2, color="gray20") +
  geom_text(aes(x=ymd("2021-03-01"), label="grocery wkrs eligible in CA", y=35), colour="gray50", angle=90, vjust = -.5, size=3)+
  geom_text(aes(x=ymd("2021-03-01")+weeks(2),  label="14 days after grocery wkrs eligible in CA", y=35), colour="gray50", angle=90, vjust = -.5, size=3)+
  scale_y_continuous("", trans="log2", limits=c(5,120))+ 
  xlab("Matched Sample")+
  scale_color_manual(values=c("cornflowerblue","red"))+
  theme_bw() +
  theme(axis.text=element_text(size=10),
        axis.title=element_text(size=10))
ggarrange(p, p1, ncol=2, nrow=1, common.legend = TRUE, legend="bottom")


```


## Figure 6: Daily changes in the Difference Between California and Oregon, Matched Sample. 

```{r fig-pt, eval=T, fig.cap="Daily changes in the Difference Between California and Oregon, $\\delta_{1t}$, Matched Sample. Error bars represent 90\\% CIs."}

matched_data1$day<-as.factor(matched_data1$Date)


PT_test_matched<-lm(ln_Newcases7~State_dum+day+State_dum:day+ln_GDP_co_2019+ln_totpop1519+perc_unemploy1519+percBA1519+urb2010+percblack1519+perclatino1519+perc_otherrace1519+perc_foreign1519+Bidenmarg +  age_65plus_prop+ indoor_dine +
                      lag3_lnNewCasescap1000+lnCumCovid1000_pre,
                     data=matched_data1)

theta<-coef(PT_test_matched)[130:244]

day=seq(from=min(ymd(matched_data1$Date))+days(1),to=max(ymd(matched_data1$Date)), by=1)
grp<-as.numeric(day>ymd("21-02-28"))
grp<-ifelse(day>ymd("21-03-14"),2,grp)
theta_means<-tapply(theta, grp, mean)

se_theta<-sqrt(diag(vcov(PT_test_matched))[130:244])

ggplot(data=cbind.data.frame(day=day,theta,se_theta,sig=as.numeric(abs(theta/se_theta)>1.645))) + geom_point(aes(x=day, y=theta))  + 
  geom_errorbar(aes(x=day, ymin=theta-1.645*se_theta, ymax=theta+1.645*se_theta), width=.1) +
  geom_hline(aes(yintercept=0))+
  geom_vline(aes(xintercept=ymd("2021-01-13")+weeks(2)), linetype=2, color="gray50") +
  geom_text(aes(x=ymd("2021-01-13")+weeks(2), label="14 days after 65+ eligible in CA", y=3), colour="gray50", angle=90, vjust = -.5, size=3)+
  geom_vline(aes(xintercept=ymd("2021-01-25")+weeks(2)), linetype=2, color="gray50") +
  geom_text(aes(x=ymd("2021-01-25")+weeks(2), label="14 days after K-12 employees eligible in OR", y=3), colour="gray50", angle=90, vjust = -.5, size=3)+
  geom_vline(aes(xintercept=ymd("2021-02-03")+weeks(2)), linetype=2, color="gray50") +
  geom_text(aes(x=ymd("2021-02-03")+weeks(2), label="14 days after the incarcerated eligible in OR", y=3), colour="gray50", angle=90, vjust = -.5, size=3)+
  geom_vline(aes(xintercept=ymd("2021-03-01")+weeks(2)), linetype=1, color="blue") +
  geom_text(aes(x=ymd("2021-03-01")+weeks(2), label="14 days after grocery wkrs eligible in CA", y=3), colour="blue", angle=90, vjust = -.5, size=3)+
  scale_y_continuous("\u03B4", limits=c(-2.5,5.5)) +
  scale_x_date("Day", limits=c(ymd("20-12-17"),ymd("21-04-11")))+
  scale_color_manual(values=c("black","red"))+
  theme_bw()+
  theme(plot.title = element_text(size = 9, face = "plain"),legend.position = "none")
```


# Appendix

## Figure 1: Centrality in the Interaction Network

```{r coleman-fig, fig.cap = "Interactions among High School Students", fig.align='center', out.width="80%", echo=TRUE, results='hide'}

data(coleman)
ggraph<-graph_from_adjacency_matrix(coleman[2,,], mode="undirected", diag=FALSE)
set.seed(1708)


ggraph <- igraph::delete.vertices(ggraph , which(igraph::degree(ggraph)==0))
degree.cent <- igraph::degree(ggraph)
degree.close<-igraph::closeness(ggraph)
degree.btw<-igraph::betweenness(ggraph)
degree.eig<-igraph::eigen_centrality(ggraph)

p<-ggplot() + geom_histogram(aes(x=degree.cent), fill="grey", colour="black")+xlab("Degree")+ylab("Count")+theme_bw()
p1<-ggplot() + geom_histogram(aes(x=degree.close), fill="grey", colour="black")+xlab("Closeness")+ylab("")+theme_bw()
p2<-ggplot() + geom_histogram(aes(x=degree.btw), fill="grey", colour="black")+xlab("Betweenness")+ylab("Count")+theme_bw()
p3<-ggplot() + geom_histogram(aes(x=degree.eig$vector), fill="grey", colour="black")+xlab("Eigenvector")+ylab("")+theme_bw()
grid.arrange(p,p1, p2,p3, nrow = 2)

```



## Figure 2: Contagion in the Interaction Network

```{r contagion-fig1, fig.cap = "Contagion in the Interaction Network", out.width="100%", fig.align='center', echo=TRUE}
par(mfrow=c(2,2), mar = c(.5, .5, 1.5, .5)) 
mynet<-(as.matrix(igraph::as_adj(ggraph)))
patient0d<-names(sort(degree.cent, decreasing=TRUE))[2]
plot_round2(mynet=mynet,patient0=patient0d,myname="Patient 0: Degree")
patient0c<-names(sort(degree.close, decreasing=TRUE))[2]
plot_round2(mynet=mynet,patient0=patient0c,myname="Patient 0: Closeness")
patient0b<-names(sort(degree.btw, decreasing=TRUE))[2]
plot_round2(mynet=mynet,patient0=patient0b,myname="Patient 0: Betweenness")
patient0e<-names(sort(degree.eig$vector, decreasing=TRUE))[2]
plot_round2(mynet=mynet,patient0=patient0e,myname="Patient 0: Eigenvector")
```


## Figure 3: Vaccination Scenarios
```{r vaccine-fig, fig.cap = "Vaccination Scenarios", out.width="100%", fig.align='center', echo=TRUE}
par(mfrow=c(2,3), mar = c(.5, .5, 1.5, .5)) 
mynet<-(as.matrix(igraph::as_adj(ggraph)))
patient0<-names(sort(degree.cent, decreasing=TRUE))[1]
#Vaccination Strategies:
n=10
vaccine_low<-as.data.frame(degree.cent) %>% filter(row.names(.)!=patient0) %>% slice_min( degree.cent,n=n, with_ties=FALSE) %>% row.names()
vaccine_centr<- as.data.frame(degree.cent) %>% filter(row.names(.)!=patient0) %>% slice_max( degree.cent,n=n, with_ties=FALSE) %>% row.names()
vaccine_close<-as.data.frame(degree.close) %>% filter(row.names(.)!=patient0) %>% slice_max( degree.close,n=n, with_ties=FALSE) %>% row.names()
vaccine_btw<-as.data.frame(degree.btw) %>% filter(row.names(.)!=patient0) %>% slice_max( degree.btw,n=n, with_ties=FALSE) %>% row.names()

degree.eig0 <- as.data.frame(degree.eig$vector) 
names(degree.eig0)<-"eig"
degree.eig0 <-degree.eig0 %>% filter(row.names(.)!=patient0)
vaccine_eig<- degree.eig0 %>% arrange(eig) %>% slice_max(eig, n=n+1, with_ties=FALSE) %>% row.names()

vulnerable<-vaccine_low


plot_round2(mynet=mynet,patient0=patient0,vaccinated=vaccine_low,vulnerable=vaccine_low,myname="Low Degree")
plot_round2(mynet=mynet,patient0=patient0,vaccinated=vaccine_centr,vulnerable=vaccine_low,myname="High Degree")
plot_round2(mynet=mynet,patient0=patient0,vaccinated=vaccine_close,vulnerable=vaccine_low,myname="High Closeness")
plot_round2(mynet=mynet,patient0=patient0,vaccinated=vaccine_btw,vulnerable=vaccine_low,myname="High Betweenness")
plot_round2(mynet=mynet,patient0=patient0,vaccinated=vaccine_eig,vulnerable=vaccine_low,myname="High Eigenvector")

```

## Table 1: ERGM Fit of the Interaction Network
```{r, eval=TRUE}


mynet<-(as.matrix(igraph::as_adj(ggraph)))
friends<-as.network(mynet, directed = FALSE)
m1<-ergm(friends~edges+ kstar(2)+gwesp(0.5, fixed = TRUE))
saveRDS(m1,"./data/ergm1.rds")


m1<-readRDS("./data/ergm1.rds") 

mytable<-cbind(coef(m1),sqrt(diag(vcov(m1))))
row.names(mytable)<-c("Edges","k-Star(2)","Gwesp (decay=0.5)")
names(mytable)<-c("Coefficient","Standard Error")

kable(mytable, caption = "ERGM Fit of the Interaction Network",digits=2) #%>% kable_styling(latex_options = "hold_position")
```

## Table 2: Round 2 Summary for 10,000 Simulations of the Interaction Network among 73 Individuals

```{r sim-res1, echo=TRUE, eval=TRUE}
source("./scripts/sum_net.R")

friends<-network(mynet, matrix.type="adjacency", 
                directed=FALSE, loops=FALSE)
m1<-readRDS("./data/ergm1.rds")
set.seed(1314)
g.sim <- simulate(m1, nsim=10000, 
                  basis=friends, control=control.simulate.ergm(
                   MCMC.burnin=1000,
                   MCMC.interval=100))

myres<-lapply(1:length(g.sim), function(i){
  sum_net(mynet=g.sim[[i]]) })
arr <- array( unlist(myres) , c(6,4,length(g.sim)) )
saveRDS(arr,"./data/sim_res1.rds")

mytable<-as.data.frame(as.matrix(NA,6,4))
for (i in 1:6){
  for (j in 1:4){
   m<-format(round(mean(arr[i,j,]),2),nsmall=2)
   se<-format(round(sd(arr[i,j,]),2),nsmall=2)
   mytable[i,j]<-paste0(m, " (",se,")")
  }
}

mytable[,2]<-c(0,10,10,10,10,10)
mytable[2,4]<-0


row.names(mytable)<-row.names(myres[[1]])
names(mytable)<-names(myres[[1]])

saveRDS(mytable,"./data/mytable1.rds")

mytable<-readRDS("./data/mytable1.rds")
mytable<-mytable[c(2,1,3,4)]
mytable %>% kable(caption = "Round 2 Summary for 10,000 Simulations of the Interaction Network among 73 Individuals",align='cccc') %>%
  kable_paper("hover", full_width = F) %>%
   kableExtra::footnote(general="Cell values are means over 10,000 simulations. Numbers in parentheses are standard deviations. \\\\textit{Not Infected} does not include \\\\textit{Vaccinated}.", escape=FALSE,threeparttable=TRUE)
```

## Figure 4: Average Infection Rates over 10,000 Simulations

```{r sim-fig, fig.cap="Average Infection Rates over 10,000 Simulations", fig.align="center"}
arr<-readRDS("./data/sim_res1.rds")

scenarios<-row.names(mytable)
groups<-names(mytable[c(2,1,3,4)])

mydata<-NULL
for (i in 1:6){
  scenario<-scenarios[i]
  for (j in c(1,3)){
   group<-groups[j]
   num<-arr[i,j,]
   d<-cbind("Scenario"=scenario,"Group"=group,num)
   mydata<-rbind(mydata,d)
  }
}


mydata %>% data.frame()  %>% mutate(num=parse_number(num),num=ifelse((Group=="Not Infected" & Scenario!="No Vaccine"),num+10, num), Group=ifelse(Group=="Not Infected","Not Infected or Vaccinated",Group)) %>% 
  group_by(Group, Scenario) %>%
  summarise(xmean=mean(num), xsd=sd(num), .groups="keep") %>%
  ggplot( ) +
    geom_point(aes(x=Scenario,y=xmean, color=Scenario),size=1)+
  geom_errorbar(aes(x=Scenario, ymin=xmean-2*xsd, ymax=xmean+2*xsd, color=Scenario),width=.5)+
  facet_grid(~Group)+theme_minimal()+theme(legend.position="bottom")+scale_color_brewer("",type="qual",palette = 3)+coord_flip()+scale_y_continuous("Count")+scale_x_discrete("Vaccine Strategy")

```

## Figure 5: Virus Spread in Networks with Cliques

```{r cliques-fig, fig.cap = "Virus Spread in Networks with Cliques", out.width="100%", fig.align="center", eval=TRUE,echo=TRUE}
set.seed(1413)
gl <- graph_from_literal(1-2,1-3,1-4,1-5,1-6,1-7,1-8,1-9,1-10,
                             2-3,2-4,2-5,2-6,2-7,2-8,2-9,2-10,
                                 3-4,3-5,3-6,3-7,3-8,3-9,3-10,
                                     4-5,4-6,4-7,4-8,4-9,4-10,
                                         5-6,5-7,5-8,5-9,5-10,
                                             6-7,6-8,6-9,6-10,
                                                 7-8,7-9,7-10,
                                                     8-9,8-10,
                                                         9-10,
                         11-12,11-13,11-14,11-15,
                               12-13,12-14,12-15,
                                     13-14,13-15,
                                           14-15,
                         10-16-11,
                         4-17-13,
                         8-17-14,
                         16-18,16-19,16-20,16-21,16-22,16-23,16-24,16-25,
                         17-26,17-27,16-28,16-29,16-30
                         )
degree.btw<-igraph::betweenness(gl)

mynet<-(as.matrix(igraph::as_adj(gl)))
patient0d<-names(sort(degree.btw, decreasing=TRUE))[1]


par(mfrow=c(2,2), mar = c(.5, .5, 1.5, .5)) 
cfg <- cluster_fast_greedy(as.undirected(gl))
cfg$membership <- c(rep(1,15), rep(2,2), rep(3,13))

V(gl)$community <-  cfg$membership
colrs <- adjustcolor( c("gray50", "green", "white"), alpha=.6)
plot(gl, vertex.color=colrs[V(gl)$community], vertex.label="")
title("Scenario 1: \n Prior to Exposure",cex.main=.8)

plot_round2(mynet=mynet,patient0=patient0d,myname="Scenario 2: \n No Vaccine")

mynet<-(as.matrix(igraph::as_adj(gl)))
patient0<-names(sort(degree.btw, decreasing=TRUE))[1]
#Vaccination Strategies:
n=5
vaccine_low<-c(1:5)
vaccine_centr<- as.data.frame(degree.btw)  %>% slice_max(degree.btw,n=n, with_ties=FALSE) %>% row.names()

vulnerable<-vaccine_low


plot_round2(mynet=mynet,patient0=patient0,vaccinated=vaccine_low,vulnerable=vaccine_low,myname="Scenario 3: Vulnerabity-Based \n Prioritization")
plot_round2(mynet=mynet,patient0=patient0,vaccinated=vaccine_centr,vulnerable=vaccine_low,myname="Scenario 4: Contact-Based \n Prioritization")


```


## Table 3: Grocery Eligibility Dates by State

```{r date-table, results='asis', eval=TRUE, echo=TRUE}
read.csv("./data/state_data1.csv", header=TRUE) %>% mutate(grocery_elig=ymd(grocery_elig)) %>% select(State, Date=grocery_elig,Gov_Rep) %>% arrange(Date) %>%  
  mutate(State=ifelse(Gov_Rep==1,paste0(State,"*"),State)) %>%
  group_by(Date) %>% 
  summarise(States=paste(State, collapse=", ")) %>%
  mutate(Date=format(Date, "%b %d"), ) %>% kable(caption = "Grocery Eligibility Dates by State",align='ll') %>%
  kable_paper("hover", full_width = F) %>%
   kableExtra::add_footnote("Republican-governed states.", notation = "symbol", escape=FALSE,threeparttable=TRUE) 
```

## Table 4: Determinants of Eligibility Date for Grocery Employees and Individuals Aged 65 and Over

```{r, eval=TRUE}
d<-read.csv("./data/state_data1.csv", header=TRUE) %>% mutate(grocery_elig=ymd(grocery_elig))
d[,22:35]<-lapply(d[,22:35], function (x) mdy(x))
d[,(ncol(d)+1):(ncol(d)+14)]<-lapply(d[,22:35], function (x) x-ymd("2020-12-01"))
d[,(ncol(d)-13):(ncol(d))]<-lapply(d[,(ncol(d)-13):(ncol(d))], function (x) as.numeric(x, "days"))

summary(m1<-coxph(Surv(num_days) ~ GDP_cap+median_age1519+log(totpop1519)+perc_unemploy1519+percBA1519+I(grocemp_pc*100)+Gov_Rep+biden_marg+urb2010+perclatino1519+percblack1519+log(Cum_covid), data= d))

summary(m2<-coxph(Surv(`Date.adults.ages.65..became.eligible.for.COVID.19.vaccination.1`) ~ GDP_cap+median_age1519+log(totpop1519)+perc_unemploy1519+percBA1519+I(grocemp_pc*100)+Gov_Rep+biden_marg+urb2010+perclatino1519+percblack1519+log(Cum_covid), data= d))

modelsummary(list("Grocery"=m1, "Age 65 Plus"=m2), estimate="{estimate}{stars} ({std.error})", statistic=NULL, 
             gof_omit='logLik|BIC|AIC',
          fmt=3,
          caption="Determinants of Eligibility Date for Grocery Employees and Individuals Aged 65 and Over",
          coef_rename=c("GDP_cap"="GDP/cap", 
                       "median_age1519"="Median Age",
                       "log(totpop1519)"="State Population, logged",
                       "perc_unemploy1519"="Unemployment Rate",
                       "percBA1519"="Percent BA Degree",
                       "I(grocemp_pc*100)"="Grocery Employees, % Population",
                       "Gov_rep"="Republican Governor",
                       "urb2010"="Urbanization",
                       "percblack1519"="Percent Black",
                       "perclatino1519"="Percent Latino",
                       "biden_marg"="Biden's Margin",
                       "log(Cum_covid)"="Cum. Covid-19 Cases, logged"))

```


## Figure 6: Decision-Making Timelines for California and Oregon. The dot–dashed horizontal line shows the Timing of the Change in the CDC Guidelines.

```{r fig-timeline, fig.cap = "Decision-Making Timelines for California and Oregon. The dot--dashed horizontal line shows the Timing of the Change in the CDC Guidelines.", fig.align="center", fig.height = 6, fig.width=5, echo=TRUE, eval=TRUE}

OR<-data.frame(State="Oregon", 
               date=c(ymd("2020-12-16"),ymd("2020-12-31"),ymd("2021-01-14"), ymd("2021-01-25"), ymd("2021-03-22")), 
               event=c("Phase 1a \n Starts","Gov. Forms \n VAC", "VAC makes \n recommendations \n for Phase 1b","Phase 1b \n Starts", "Grocery \n Workers  \n Eligible"), 
               disloc=c(.5,.5,4,1,1))

CA<-data.frame(State="California", 
               date=c(ymd("2020-11-23"), ymd("2020-12-16"), ymd("2020-12-23"),ymd("2021-01-04"),ymd("2021-01-13"),ymd("2021-03-01")), 
               event=c("Gov. Forms \n VAC","Phase 1a \n Starts","VAC makes \n recommendations \n for Phase 1b","Gov. Announces \n Make-up of Phases 1b and 1c", "Phase 1b \n Starts", "Grocery \n Workers \n Eligible"),
               disloc=c(-1,-.5,-3,-5,-.5,-.5) )

data <- rbind.data.frame(OR,CA)  %>% arrange(date) %>% mutate(X=as.numeric(ymd("2020-11-23")- date , "days"), 
                                                              disloc1=ifelse(State=="Oregon", disloc+3.5, disloc))


myb<-sort(c(unique(data$X)[-c(6:7)],-50), decreasing =T)
myl<-sort(c(unique(data$date)[-c(6:7)],ymd("2021-01-12"))) 

ggplot() +
    geom_segment(data=data, aes(x = X,y = disloc, xend = X, yend = 0, color=State)) +
    geom_segment(data=data, aes(x =3, y = 0, xend = min(X)-10, yend = 0), arrow = arrow(length = unit(x = 0.2,units = 'cm'),type = 'closed')) +
    geom_text(data=data, aes(x = X,y = disloc1,label = event, color=State), show.legend=FALSE, size=3.5,lineheight = .8,hjust = 1.0,vjust = .8,parse = FALSE) +
    geom_point(data=data, aes(x = X,y = disloc, color=State)) +
    geom_segment(aes(x = -50,y = -8, xend = -50,yend = 7), linetype=6, color="black", size=.5) +
   geom_text(aes(x=-57, y=-8.5, label="Change \n in CDC \n Guidelines"), size=3.5, color="black") +
    scale_x_continuous("", breaks = myb,labels = format(myl,format="%b-%d")) +
    scale_y_continuous("", limits=c(-9,8)) +
    theme_bw() +
    coord_flip()+
    theme(axis.text.x = element_blank(),axis.text.y = element_text(size=10),axis.ticks = element_blank(),axis.title.x = element_blank(),axis.title.y = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.border = element_blank())+
  guides(fill = guide_legend(reverse = TRUE))+
  theme(legend.position='bottom', legend.direction='horizontal') 

```

## Figure 7: Prioritization Timelines for California and Oregon. The dashed vertical line shows the data of open vaccine eligibility in both states.

```{r fig-priorities, fig.cap = "Prioritization Timelines for California and Oregon. The dashed vertical line shows the data of open vaccine eligibility in both states.", fig.align="center", out.height="100%", echo=TRUE, eval=TRUE}
d<-read.csv("./data/state_data1.csv", header=TRUE)
d[,22:35]<-lapply(d[,22:35], function (x) mdy(x))
d[,(ncol(d)+1):(ncol(d)+14)]<-lapply(d[,22:35], function (x) x-ymd("2020-12-01"))
d[,(ncol(d)-13):(ncol(d))]<-lapply(d[,(ncol(d)-13):(ncol(d))], function (x) as.numeric(x, "days"))

mynames<-names(d)
mynames[22:35]<-paste0("Date ",c("Age 80 Plus","Age 75 Plus","Age 70 Plus","Age 65 Plus","Age 60 Plus","Age 55 Plus","Age 50 Plus","Age 45 Plus","Age 40 Plus","Age 30 Plus","K-12","Grocery and Agriculture","Incarcerated","General Public"))
names(d)<-mynames
d[,(ncol(d)+1):(ncol(d)+14)]<-lapply(d[,22:35], function (x) x- ymd("2021-04-19"))
d[,(ncol(d)-13):(ncol(d))]<-lapply(d[,(ncol(d)-13):(ncol(d))], function (x) as.numeric(x, "days"))

mynames<-names(d)
mynames[(ncol(d)-13):(ncol(d))]<-c("Age 80 Plus","Age 75 Plus","Age 70 Plus","Age 65 Plus","Age 60 Plus","Age 55 Plus","Age 50 Plus","Age 45 Plus","Age 40 Plus","Age 30 Plus","K-12","Grocery and Agriculture","Incarcerated","General Public")
names(d)<-mynames


mylabs<-subset(d, State %in% c("California","Oregon")) %>% select(State,"Date Age 80 Plus", "Date Age 75 Plus", "Date Age 70 Plus","Date Age 65 Plus","Date Age 60 Plus","Date Age 55 Plus","Date Age 50 Plus","Date Age 45 Plus","Date Age 40 Plus","Date Age 30 Plus","Date K-12","Date Grocery and Agriculture", "Date Incarcerated","Date General Public")

subset(d, State %in% c("California","Oregon")) %>% select(State,"General Public","Age 50 Plus","Grocery and Agriculture","Incarcerated","K-12","Age 65 Plus","Age 70 Plus","Age 75 Plus","Age 80 Plus") %>% 
  pivot_longer(cols=2:10, values_to="Days", names_to="Group") %>%
  mutate(Group=factor(Group, levels = c("General Public","Age 50 Plus","Incarcerated","Grocery and Agriculture","K-12","Age 65 Plus","Age 70 Plus","Age 75 Plus","Age 80 Plus")), State=factor(State, levels=c("Oregon", "California"))) %>%
  ggplot()+ geom_bar(aes(x=Group, y=Days-30, fill=State), stat="identity", position = "dodge")+
  geom_hline(yintercept=-30, linetype=2, size=1, color="black")+ 
  scale_y_continuous("",limits=c(-145,0),breaks=c(-138,-107,-79,-48,-18),labels=format(c(ymd("2021-01-01"),ymd("2021-02-01"),ymd("2021-03-01"),ymd("2021-04-01"),ymd("2021-05-01")), format="%b-%d"))+
  scale_x_discrete("")+
  scale_fill_manual(values=c("#00BFC4","#F8766D"))+
  coord_flip(clip = "off")+
  theme_bw()+
  guides(fill = guide_legend(reverse = TRUE))+
  theme(legend.position='bottom', legend.direction='horizontal') 
```

##Table 5: The Effect of Vaccine Eligibility to Grocery Employees on Covid-19 Hospitalizations (logged)

```{r hosp, eval=TRUE}
mydata<-read.csv("./data/CAOR_data.csv", header=TRUE) %>% mutate(Date=ymd(Date)) %>% filter(Date<ymd("2021-03-29")+weeks(4) & Date>ymd("2020-12-16")) %>% select(-c("grocery_day"))  

mydata<- mydata %>% filter(Date==ymd("2020-12-17")) %>% 
  select(State, County, lnCumHosp1000_pre=lnCumHosp1000,CumHosp_pre=CumHosp,CumCovid_pre=Cum_covid,lnCumCovid1000_pre=lnCumCovid1000) %>%  
  left_join(mydata, by=c("State","County"))


tomatch<-mydata %>% filter(Date==ymd("2020-12-17")) %>% select(State, County,GDP_co_2019,ln_GDP_co_2019,ln_totpop1519,totpop1519,perc_unemploy1519,percBA1519,urb2010,percblack1519,perc_otherrace1519,perclatino1519,perc_foreign1519,Bidenmarg, age_65plus_prop, indoor_dine,
                           CumCovid_pre,lnCumCovid1000_pre) 

todrop<-c("State","County","totpop1519","GDP_co_2019","CumCovid_pre","grocery_wkrs100", "age_65plus_prop", "indoor_dine","CumHosp_pre","lnCumHosp1000_pre")
mat1 <- cem(treatment = "State",
            cutpoints=list(ln_GDP_co_2019=4,
                           ln_totpop1519=6,
                           perc_unemploy1519=2,
                           percBA1519=2, 
                           urb2010=2,
                           percblack1519=3,
                           perc_otherrace1519=3,
                           perclatino1519=2, 
                           perc_foreign1519=2,
                          Bidenmarg=4),
            baseline.group = "Oregon",
            drop=todrop,data = tomatch, keep.all = TRUE)

matches<-tomatch[mat1$matched==TRUE,] %>% mutate(matched="TRUE") %>% select(State,County,matched)
matched_data1<-mydata %>% left_join(matches, by=c("State","County")) %>% filter(matched==TRUE)

mydata <- mydata %>% mutate(grocery=ifelse(State=="California" & Date>=ymd("2021-03-01")+weeks(4), 1,0),
         grocery_day=time_length(interval(ymd("2021-03-01")+weeks(4),Date),unit="day"),
         grocery_day=ifelse(Date<ymd("2021-03-01")+weeks(4),0,grocery_day),
         State_grocery_day=State_dum*grocery_day)

matched_data1 <- matched_data1 %>% mutate(grocery=ifelse(State=="California" & Date>=ymd("2021-03-01")+weeks(4), 1,0),
         grocery_day=time_length(interval(ymd("2021-03-01")+weeks(4),Date),unit="day"),
         grocery_day=ifelse(Date<ymd("2021-03-01")+weeks(4),0,grocery_day),
         State_grocery_day=State_dum*grocery_day)


m1_all_hosp<-lm(ln_Hosp7~grocery_day+State_dum+State_grocery_day+lnCumHosp1000_pre+ln_GDP_co_2019+ln_totpop1519+perc_unemploy1519+percBA1519+urb2010+percblack1519+perclatino1519+perc_otherrace1519+perc_foreign1519+Bidenmarg+ age_65plus_prop+ indoor_dine+lag14_lnNewCasescap1000, data=mydata)

m1_matched_hosp<-lm(ln_Hosp7~grocery_day+State_dum+State_grocery_day+lnCumHosp1000_pre+ln_GDP_co_2019+ln_totpop1519+perc_unemploy1519+percBA1519+urb2010+percblack1519+perclatino1519+perc_otherrace1519+perc_foreign1519+Bidenmarg+ age_65plus_prop+ indoor_dine+lag14_lnNewCasescap1000, data=matched_data1)


modelsummary(list("Full Sample"=m1_all_hosp, "Matched Sample"=m1_matched_hosp), estimate="{estimate}{stars} ({std.error})", statistic=NULL, 
             gof_omit='logLik|BIC|AIC',
          fmt=3,
          caption="The Effect of Vaccine Eligibility to Grocery Employees on Covid-19 Hospitalizations (logged)",
          coef_rename=c("grocery_day"="Day of Treatment", 
                       "State_dum"="California",
                       "State_grocery_day"="California*Day of Treatment",
                       "lnCumHosp1000_pre"="Cumulative Hospitalizations, logged",
                       "ln_GDP_co_2019"="County GDP, logged",
                       "ln_totpop1519"="County Population, logged",
                       "perc_unemploy1519"="Unemployment Rate",
                       "percBA1519"="Percent BA Degree",
                       "urb2010"="Urbanization",
                       "percblack1519"="Percent Black",
                       "perclatino1519"="Percent Latino",
                       "perc_otherrace1519"="Percent Other Race",
                       "perc_foreign1519"="Percent Foreign",
                       "Bidenmarg"="Biden's Margin",
                       "grocery_wkrs100"="Grocery Wkrs, per 100 res.",
                       "age_65plus_prop"="Prop. Aged 65+",
                       "indoor_dine"="Indoor Dining Ban",
                       "lag3_lnHosp1000"="Lagged DV",
                       "(Intercept)"="Constant"))

```


## Figure 8: Marginal Effect of Vaccine Eligibility to Grocery Workers On Hospitalizations. Error bars represent 95% CIs.

```{r fig-marg-hosp, fig.cap="Marginal Effect of Vaccine Eligibility to Grocery Workers On Hospitalizations. Error bars represent 95\\% CIs.", eval=TRUE, echo=TRUE}

grocery_day=seq(from=min(mydata$grocery_day),to=max(mydata$grocery_day), by=1)

marg_state<-m1_all_hosp$coef[3]+m1_all_hosp$coef[4]*grocery_day
se_state<-sqrt(vcov(m1_all_hosp)[3,3]+vcov(m1_all_hosp)[4,4]*(grocery_day^2)+2*grocery_day*vcov(m1_all_hosp)[3,4])

p<-ggplot(data=cbind.data.frame(grocery_day,marg_state,se_state)) + geom_point(aes(x=grocery_day, y=marg_state))  + 
  geom_errorbar(aes(x=grocery_day, ymin=marg_state-2*se_state, ymax=marg_state+2*se_state), width=.1) +
  geom_hline(aes(yintercept=0))+
 scale_y_continuous("Marginal Effect of California", limits=c(-.7,.65), breaks=seq(from=-.7,to=.7, by=.2)) +
  ggtitle("Full Sample (Model 1)")+
  scale_x_continuous("Day of Treatment", limits=c(0,20))+
  theme_bw()+
  theme(plot.title = element_text(size = 9, face = "plain"))

# (Matched Sample):
grocery_day=seq(from=min(matched_data1$grocery_day),to=max(matched_data1$grocery_day), by=1)

marg_state<-m1_matched_hosp$coef[3]+m1_matched_hosp$coef[4]*grocery_day
se_state<-sqrt(vcov(m1_matched_hosp)[3,3]+vcov(m1_matched_hosp)[4,4]*(grocery_day^2)+2*grocery_day*vcov(m1_matched_hosp)[3,4])

p1<-ggplot(data=cbind.data.frame(grocery_day,marg_state,se_state)) + geom_point(aes(x=grocery_day, y=marg_state))  + 
  geom_errorbar(aes(x=grocery_day, ymin=marg_state-2*se_state, ymax=marg_state+2*se_state), width=.1) +
  geom_hline(aes(yintercept=0))+
  scale_y_continuous("", limits=c(-.7,.65), breaks=seq(from=-.7,to=.7, by=.2)) +
  ggtitle("Matched Sample (Model 2)")+
  scale_x_continuous("Day of Treatment", limits=c(0,20))+
  theme_bw()+
  theme(plot.title = element_text(size = 9, face = "plain"))

myfig<-ggarrange(p+ rremove("ylab") + rremove("xlab"),
                 p1+ rremove("ylab") + rremove("xlab"),
                 ncol=2, nrow=1, common.legend = TRUE, legend="bottom")
annotate_figure(myfig, left = "Marginal Effect of California",
                bottom = "Day of Treatment")

```


## Table 6: The Effect of Vaccine Eligibility to Grocery Employees on Covid-19 Deaths (logged)

```{r deaths, eval=TRUE}

mydata<-read.csv("./data/CAOR_data.csv", header=TRUE) %>% mutate(Date=ymd(Date)) %>% filter(Date>ymd("2020-12-16"))
mydata<- mydata %>% filter(Date==ymd("2020-12-17")) %>% 
  select(State, County, lnCumDeaths1000_pre=lnCumDeaths1000,Cum_deaths_pre=Cum_deaths,CumCovid_pre=Cum_covid,lnCumCovid1000_pre=lnCumCovid1000) %>%  
  left_join(mydata, by=c("State","County"))

tomatch<-mydata %>% filter(Date==ymd("2020-12-17")) %>% select(State, County,GDP_co_2019,ln_GDP_co_2019,ln_totpop1519,totpop1519,perc_unemploy1519,percBA1519,urb2010,percblack1519,perc_otherrace1519,perclatino1519,perc_foreign1519,Bidenmarg, age_65plus_prop, indoor_dine,
                           CumCovid_pre,lnCumCovid1000_pre) 

todrop<-c("State","County","totpop1519","GDP_co_2019","CumCovid_pre","grocery_wkrs100", "age_65plus_prop", "indoor_dine","CumHosp_pre","lnCumHosp1000_pre")
mat1 <- cem(treatment = "State",
            cutpoints=list(ln_GDP_co_2019=4,
                           ln_totpop1519=6,
                           perc_unemploy1519=2,
                           percBA1519=2, 
                           urb2010=2,
                           percblack1519=3,
                           perc_otherrace1519=3,
                           perclatino1519=2, 
                           perc_foreign1519=2,
                          Bidenmarg=4),
            baseline.group = "Oregon",
            drop=todrop,data = tomatch, keep.all = TRUE)




matches<-tomatch[mat1$matched==TRUE,] %>% mutate(matched="TRUE") %>% select(State,County,matched)
matched_data1<-mydata %>% left_join(matches, by=c("State","County")) %>% filter(matched==TRUE)


mydata <- mydata %>% mutate(grocery=ifelse(State=="California" & Date>=ymd("2021-03-01")+weeks(6),1,0),    grocery_day=time_length(interval(ymd("2021-03-01")+weeks(6),Date),unit="day"),
         grocery_day=ifelse(Date<ymd("2021-03-01")+weeks(6),0,grocery_day), 
         State_grocery_day=State_dum*grocery_day)

matched_data1 <- matched_data1 %>% mutate(grocery=ifelse(State=="California" & Date>=ymd("2021-03-01")+weeks(6),1,0),    grocery_day=time_length(interval(ymd("2021-03-01")+weeks(6),Date),unit="day"),
         grocery_day=ifelse(Date<ymd("2021-03-01")+weeks(6),0,grocery_day), 
         State_grocery_day=State_dum*grocery_day)


m1_all_deaths<-lm(ln_Deaths7~grocery_day+State_dum+State_grocery_day+lnCumDeaths1000_pre+ln_GDP_co_2019+ln_totpop1519+perc_unemploy1519+percBA1519+urb2010+percblack1519+perclatino1519+perc_otherrace1519+perc_foreign1519+Bidenmarg+ age_65plus_prop+ indoor_dine
                  +lag28_lnNewCasescap1000, 
                  data=mydata)

m1_matched_deaths<-lm(ln_Deaths7~grocery_day+State_dum+State_grocery_day
                      +lnCumDeaths1000_pre+ln_GDP_co_2019+ln_totpop1519+perc_unemploy1519+percBA1519+urb2010+percblack1519+perclatino1519+perc_otherrace1519+perc_foreign1519+Bidenmarg+ age_65plus_prop+ indoor_dine+lag28_lnNewCasescap1000
                      , data=matched_data1)


modelsummary(list("Full Sample"=m1_all_deaths, "Matched Sample"=m1_matched_deaths), estimate="{estimate}{stars} ({std.error})", statistic=NULL, 
             gof_omit='logLik|BIC|AIC',
          fmt=3,
          caption="The Effect of Vaccine Eligibility to Grocery Employees on Covid-19 Deaths (logged)",
          coef_rename=c("grocery_day"="Day of Treatment", 
                       "State_dum"="California",
                       "State_grocery_day"="California*Day of Treatment",
                       "lnCumDeaths1000_pre"="Cumulative Deaths, logged",
                       "ln_GDP_co_2019"="County GDP, logged",
                       "ln_totpop1519"="County Population, logged",
                       "perc_unemploy1519"="Unemployment Rate",
                       "percBA1519"="Percent BA Degree",
                       "urb2010"="Urbanization",
                       "percblack1519"="Percent Black",
                       "perclatino1519"="Percent Latino",
                       "perc_otherrace1519"="Percent Other Race",
                       "perc_foreign1519"="Percent Foreign",
                       "Bidenmarg"="Biden's Margin",
                       "grocery_wkrs100"="Grocery Wkrs, per 100 res.",
                       "age_65plus_prop"="Prop. Aged 65+",
                       "indoor_dine"="Indoor Dining Ban",
                       "lag3_lnHosp1000"="Lagged DV",
                       "(Intercept)"="Constant"))

```

## Figure 9: Marginal Effect of Vaccine Eligibility to Grocery Workers On Covid-19 Deaths. Error bars represent 95% CIs.


```{r fig-marg-death, fig.cap="Marginal Effect of Vaccine Eligibility to Grocery Workers On Covid-19 Deaths. Error bars represent 95\\% CIs.", eval=TRUE}
grocery_day=seq(from=min(mydata$grocery_day),to=max(mydata$grocery_day), by=1)

marg_state<-m1_all_deaths$coef[3]+m1_all_deaths$coef[4]*grocery_day
se_state<-sqrt(vcov(m1_all_deaths)[3,3]+vcov(m1_all_deaths)[4,4]*(grocery_day^2)+2*grocery_day*vcov(m1_all_deaths)[3,4])

p<-ggplot(data=cbind.data.frame(grocery_day,marg_state,se_state)) + geom_point(aes(x=grocery_day, y=marg_state))  + 
  geom_errorbar(aes(x=grocery_day, ymin=marg_state-2*se_state, ymax=marg_state+2*se_state), width=.1) +
  geom_hline(aes(yintercept=0))+
  scale_y_continuous("Marginal Effect of California", limits=c(-.17,.25)) +
  ggtitle("Full Sample (Model 1)")+
  scale_x_continuous("Day of Treatment", limits=c(0,20))+
  theme_bw()+
  theme(plot.title = element_text(size = 9, face = "plain"))


grocery_day=seq(from=min(matched_data1$grocery_day),to=max(matched_data1$grocery_day), by=1)

marg_state<-m1_matched_deaths$coef[3]+m1_matched_deaths$coef[4]*grocery_day
se_state<-sqrt(vcov(m1_matched_deaths)[3,3]+vcov(m1_matched_deaths)[4,4]*(grocery_day^2)+2*grocery_day*vcov(m1_matched_deaths)[3,4])

p1<-ggplot(data=cbind.data.frame(grocery_day,marg_state,se_state)) + geom_point(aes(x=grocery_day, y=marg_state))  + 
  geom_errorbar(aes(x=grocery_day, ymin=marg_state-2*se_state, ymax=marg_state+2*se_state), width=.1) +
  geom_hline(aes(yintercept=0))+
  scale_y_continuous("", limits=c(-.17,.25)) +
  ggtitle("Matched Sample (Model 2)")+
  scale_x_continuous("Day of Treatment", limits=c(0,20))+
  theme_bw()+
  theme(plot.title = element_text(size = 9, face = "plain"))

myfig<-ggarrange(p+ rremove("ylab") + rremove("xlab"),
                 p1+ rremove("ylab") + rremove("xlab"),
                 ncol=2, nrow=1, common.legend = TRUE, legend="bottom")
annotate_figure(myfig, left = "Marginal Effect of California",
                bottom = "Day of Treatment")

```

## Table 7: The Effect of Vaccine Eligibility to Grocery Employees on New Daily Covid-19 Cases (logged)

```{r}
mydata<-read.csv("./data/CAOR_data.csv", header=TRUE) %>% mutate(Date=ymd(Date)) %>% filter(Date<ymd("2021-03-29")+weeks(2) & Date>ymd("2020-12-16"))

mydata %>% filter(Date==ymd("2020-12-17")) %>% 
  select(State, County, lnCumCovid1000_pre=lnCumCovid1000,CumCovid_pre=Cum_covid) %>%  
  left_join(mydata, by=c("State","County")) -> mydata


tomatch<-mydata %>%  filter(Date==ymd("2020-12-17")) %>% select(State, County, GDP_co_2019,ln_GDP_co_2019,ln_totpop1519,totpop1519,perc_unemploy1519,percBA1519,urb2010,percblack1519,perc_otherrace1519,perclatino1519,perc_foreign1519,Bidenmarg, grocery_wkrs100, age_65plus_prop, indoor_dine,
                           CumCovid_pre,lnCumCovid1000_pre) 
todrop<-c("State","County","totpop1519","GDP_co_2019","CumCovid_pre","grocery_wkrs100", "age_65plus_prop", "indoor_dine")
mat1 <- cem(treatment = "State",
            cutpoints=list(ln_GDP_co_2019=4,
                           ln_totpop1519=6,
                           perc_unemploy1519=2,
                           percBA1519=2, 
                           urb2010=2,
                           percblack1519=3,
                           perc_otherrace1519=3,
                           perclatino1519=2, 
                           perc_foreign1519=2,
                          Bidenmarg=4),
            baseline.group = "Oregon",
            drop=todrop,data = tomatch, keep.all = TRUE)

matches<-tomatch[mat1$matched==TRUE,] %>% mutate(matched="TRUE") %>% select(State,County,matched)
matched_data1<-mydata %>% left_join(matches, by=c("State","County")) %>% filter(matched==TRUE)

mydata$State_grocery_day<-mydata$State_dum*mydata$grocery_day
matched_data1$State_grocery_day<-matched_data1$State_dum*matched_data1$grocery_day
m1_all_fe<-lm(ln_Newcases7~grocery_day+State_dum+State_grocery_day+lag3_lnNewCasescap1000+factor(County), data=mydata)

m1_matched_fe<-lm(ln_Newcases7~grocery_day+State_dum+State_grocery_day+lag3_lnNewCasescap1000+factor(County), data=matched_data1)

modelsummary(list("Full Sample"=m1_all_fe, "Matched Sample"=m1_matched_fe), estimate="{estimate}{stars} ({std.error})", statistic=NULL, 
             gof_omit='logLik|BIC|AIC',
          fmt=3,
          caption="The Effect of Vaccine Eligibility to Grocery Employees on New Daily Covid-19 Cases (logged)",
          coef_rename=c("grocery_day"="Day of Treatment", 
                       "State_dum"="California",
                       "State_grocery_day"="California*Day of Treatment",
                       "(Intercept)"="Constant"))
```

